function testspec()
%
   h = findobj('Tag', 'testspec');
   if(isempty(h))
      h = figure('Tag', 'testspec', 'name', 'test spectral analysis', 'NumberTitle', 'off');
   end;
   figure(h);
   clf;
   
fs = 10000; % sample at 10 kHz
RATES = 1/fs; % rate
flomb = [1:500]; % express F in Hz
X=[0:RATES:0.1]; % and time in sec
for i = 1:length(X)
   N(i) = randn;
end;
Y=sin(2*pi*X/0.01)+0*N; % corresponds to 10 msec period = 100 Hz; can add noise too.
subplot('Position', [0.1, 0.55, 0.8, 0.4]);
plot(X,Y);

[Pv, f] = pwelch(Y, 1024, fs, 128, 32);
Pv = Pv(1:end,1)*fs;  % note unscaling by Fs here... to undo what pwelch does...
[fl, Pvl, conf] = lombpsd(X, Y, f, 0.05); % compute lomb and get confidence level
Pv = Pv/((2 * var(Y))); % rescale fft version in same way as lomb
subplot('Position', [0.1, 0.1, 0.8, 0.4]);
plot (log10(f(2:end)), Pv(2:end), '-b');
hold on;
plot(log10(fl(2:end)), Pvl(2:end), '-r');
plot(log10([fl(2) fl(end)]), [conf conf], '--k');
